
import pandas as pd
from scipy.stats import entropy
import numpy as np

# For visualisation
import seaborn as sns
import matplotlib.pyplot as plt

# For Anomaly Detection
from sklearn.ensemble import IsolationForest

# For Oracle Correlation
from scipy.stats import pearsonr
import random

# For Rasterplot
from rastermap import Rastermap
from scipy.stats import zscore

# For AllenSDK
import os
from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache

# Filter Ignore
# import warnings
# warnings.filterwarnings('ignore')
from collections import OrderedDict, defaultdict

import numpy as np
import torch
import torch.nn as nn
from itertools import product
from torch.nn import functional as F
import numpy as np
import scipy.signal
import torch
import math
from math import ceil
from torch import nn as nn
# from .module import Module
from torch.nn import Parameter
from torch.nn import functional as F
from torch.nn.init import xavier_normal
from torch.nn.modules.utils import _pair
from torch.utils.data import DataLoader
from tqdm import tqdm
import matplotlib.pyplot as plt
import scipy
import pandas as pd
import torchvision
from torchvision.transforms import Resize

import os
import pandas as pd
from PIL import Image
import torch
from torch.utils.data import Dataset

from torchvision.transforms import Resize
from numpy.fft import fft2, fftshift, fftfreq


data_dir = "./allendata"
#Again, we might want to use another folder if we do want to overwrite potential data we have
#data_dir = "./converted"

#From here on everything works with both AllenSDK and the MiniSDK!
manifest_path = os.path.join(data_dir, "manifest.json")
cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)
sessions = cache.get_session_table() #Returns a Pandas dataframe
sessions.head() #Sneak peek of the table
#Grab our (filtered) data from our favorite session. We take a female mice with nice unit count
session_id = 798911424
oursession = cache.get_session_data(session_id, timeout=3000)


image_id = oursession.get_stimulus_table("natural_scenes").index.values
visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values

#Tables of the spikes happened during these presentations
visam_spikes = oursession.presentationwise_spike_times(
    stimulus_presentation_ids=image_id,
    unit_ids = visam_ids
)
visam_spikes


nat_scenes_temp = oursession.get_stimulus_table('natural_scenes')
#nat_scenes_temp.groupby('frame').count()
nat_scenes_temp


nat_scenes_temp = nat_scenes_temp.reset_index()
df = nat_scenes_temp.copy()
# Assuming df is your DataFrame with columns: 'stimuli', 'units', 'time'
# Sort the DataFrame by 'time' column
df_sorted = df.sort_values(by='start_time')
df_sorted.reset_index(inplace = True)
#final_data = df_first_presentation.merge(visam_spikes, on = 'stimulus_presentation_id', how = 'inner')
final_data = df_sorted.merge(visam_spikes, on = 'stimulus_presentation_id', how = 'inner')
final_data


#Get the ids of the images and the units
#image_id = oursession.get_stimulus_table("natural_scenes").index.values
#visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values
image_id = final_data['stimulus_presentation_id'].values
visam_ids = final_data['unit_id'].values
print('done')


#Get the spikes during the session
#visam_spikes = oursession.presentationwise_spike_times(
#    stimulus_presentation_ids=image_id,
#    unit_ids = visam_ids
#)

visam_spikes = visam_spikes[visam_spikes['stimulus_presentation_id'].isin(image_id)]
visam_spikes = visam_spikes[visam_spikes['unit_id'].isin(visam_ids)]
print('done2')

#Add a new column to our table and fill it by counting number of rows with a stimulus presentation and unitid
visam_spikes["count"] = np.zeros(len(visam_spikes))
visam_spikes = visam_spikes.groupby(["stimulus_presentation_id", "unit_id"]).count()

#Employ "pivot table" to use the information from visam_spikes to generate a new table where
#we indicate index, columns, and fill it with the count.
spikes_presentation = pd.pivot_table(
    visam_spikes,
    values="time_since_stimulus_presentation_onset",
    index="stimulus_presentation_id",
    columns="unit_id",
    fill_value=0.0,
    aggfunc=np.sum
)

#This is 50 times faster!
spikes_presentation

done
done2


titles = final_data[['stimulus_presentation_id', 'frame']]
final_cnn = spikes_presentation.merge(titles, on = 'stimulus_presentation_id').drop_duplicates()
final_cnn.head()


final_cnn = final_cnn.groupby('frame').mean()
final_cnn


final_cnn.to_csv('dataset_for_agne_mean.csv')


df = pd.read_csv("dataset_for_agne_mean.csv")
df


data_dir = "./allendata"
manifest_path = os.path.join(data_dir, "manifest.json")
cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)


sessions = cache.get_session_table() #Returns a Pandas dataframe
session_id = 798911424
oursession = cache.get_session_data(session_id, timeout=3000)


nat_scenes_temp = oursession.get_stimulus_table('natural_scenes')
nat_scenes_temp = nat_scenes_temp.reset_index()
nat_scenes_temp


image_id = oursession.get_stimulus_table("natural_scenes").index.values
visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values

visam_spikes = oursession.presentationwise_spike_times(
    stimulus_presentation_ids=image_id,
    unit_ids = visam_ids
)
visam_spikes.reset_index(inplace=True)
visam_spikes


# We're interested in the spikes that are in response to the natural scenes
df_new = visam_spikes.copy()
# Drop the columns that we don't need
df_new.drop(columns=["time_since_stimulus_presentation_onset"], inplace=True)
# Count the number of spikes per stimulus presentation
df_new = df_new.groupby(['stimulus_presentation_id', 'unit_id']).count()
df_new.reset_index(inplace=True)
df_new


# Take natural scenes in order to take also the image IDs
df_natural = oursession.get_stimulus_table('natural_scenes')
df_natural.reset_index(inplace = True)
# Drop useless columns
df_natural.drop(columns=['stimulus_block', 'start_time', 'stop_time', 'stimulus_name', 'duration', 'stimulus_condition_id'], inplace=True)
# Rename the columns
df_natural.rename(columns={'stimulus_presentation_id': 'stimulus_presentation_id', 'frame':'image_id'}, inplace=True)
df_natural


# Merge the two datasets
final_df = df_new.merge(df_natural, how ='left', on='stimulus_presentation_id' )
final_df


# Rename the columns
final_df.rename(columns={'stimulus_presentation_id': "trial_id", "unit_id": "neuron_id", "spike_time": "spike_count", "image_id": "image_id"}, inplace = True)
# Convert image_id to integer
convert_dict = {'image_id':int}
final_df = final_df.astype(convert_dict)


final_df


final_df.to_csv("oracle_df.csv", index=False)


min_time = visam_spikes['spike_time'].min()
max_time = visam_spikes['spike_time'].max()
binsize = 0.1 # BINSIZE (!!!!)
bins = np.arange(min_time, max_time, binsize)
min_time, max_time

(5909.794764092189, 8569.517004859086)


def count_spikes(group):
    counts, _ = np.histogram(group['spike_time'], bins=bins)
    return counts

spike_counts = visam_spikes.groupby('unit_id').apply(count_spikes)

# Convert the Series of arrays into a 2D NumPy array
spike_counts_matrix = np.stack(spike_counts.values)
spike_counts_matrix

array([[2, 3, 2, ..., 2, 1, 2],
       [1, 1, 2, ..., 0, 0, 1],
       [1, 1, 3, ..., 2, 2, 2],
       ...,
       [1, 1, 2, ..., 1, 3, 3],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 2, ..., 0, 0, 0]])


spike_counts_zscored = zscore(spike_counts_matrix, axis=1)
spike_counts_zscored.shape

(135, 26597)


model = Rastermap(n_clusters=None, # None turns off clustering and sorts single neurons
                  n_PCs=64, # use fewer PCs than neurons
                  locality=0.15, # some locality in sorting (this is a value from 0-1)
                  time_lag_window=15, # use future timepoints to compute correlation
                  grid_upsample=0, # 0 turns off upsampling since we're using single neurons
                ).fit(spike_counts_zscored)
y = model.embedding # neurons x 1
isort = model.isort

normalizing data across axis=1
projecting out mean along axis=0
data normalized, 0.06sec
sorting activity: 135 valid samples by 26597 timepoints
n_PCs = 64 computed, 4.75sec
skipping clustering, n_clusters is None
clusters sorted, time 135.09sec
rastermap complete, time 135.10sec


isort = model.isort
embedding = model.embedding


xmin = 22000

xmax = 28000

fig, ax = plt.subplots(figsize=(12, 8), dpi=200)

im = ax.imshow(spike_counts_zscored[isort, xmin:xmax], cmap="gray_r", vmin=0, vmax=1.2, aspect="auto")
ax.set_xlabel("Time (bins)")
ax.set_ylabel("Neurons (sorted)")
ax.set_title("Raster Plot of Neural Activity")

cbar = plt.colorbar(im, ax=ax, pad=0.01)
cbar.set_label('Spike count')

plt.show()


# Define different xmin and xmax for each plot
xmin = [0, 8000, 22000]
xmax = [4800, 12000, 28000]

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 6), dpi=200)

for i, ax in enumerate(axes.flat):
    # Plot each segment of the data
    im = ax.imshow(spike_counts_zscored[isort, xmin[i]:xmax[i]], cmap="gray_r", vmin=0, vmax=1.2, aspect="auto")
    ax.set_xlabel("Time (bins)")
    ax.set_ylabel("Neurons (sorted)")
    ax.set_title(f"Raster Plot from {xmin[i]} to {xmax[i]} bins")


plt.tight_layout()
plt.show()


data = pd.read_csv('dataset_for_agne_mean.csv')
data.drop(['frame', 'stimulus_presentation_id'], axis=1, inplace=True)
data = data.transpose()


correlation_matrix = data.transpose().corr()
# Plotting the correlation matrix as a heatmap
plt.figure(figsize=(10, 10))  # Adjust the size as needed
sns.heatmap(correlation_matrix, cmap='coolwarm', center=0)
plt.show()


np.fill_diagonal(correlation_matrix.values, np.nan)


highest_correlations = correlation_matrix.max()
lowest_correlation = correlation_matrix.min()
most_correlated_neurons = correlation_matrix.idxmax()
least_correlated_neurons = correlation_matrix.idxmin()
results = pd.DataFrame({
    'Neuron': highest_correlations.index,
    'Highest Correlation': highest_correlations.values,
    'Most Correlated Neuron': most_correlated_neurons.values,
    'Lowest Correlation': lowest_correlation.values,
    'Least Correlated Neuron': least_correlated_neurons.values
})
print(results)

        Neuron  Highest Correlation Most Correlated Neuron  \
0    951092050             0.297214              951093080   
1    951092075             0.309912              951098279   
2    951092303             0.416028              951093160   
3    951092369             0.256688              951098530   
4    951092398             0.655551              951092410   
..         ...                  ...                    ...   
130  951099272             0.342826              951098414   
131  951099306             0.797598              951094081   
132  951099320             0.462756              951098242   
133  951099475             0.392041              951098430   
134  951099491             0.403545              951092303   

     Lowest Correlation Least Correlated Neuron  
0             -0.242279               951092488  
1             -0.239617               951093510  
2             -0.476467               951092949  
3             -0.240787               951093608  
4             -0.418246               951092303  
..                  ...                     ...  
130           -0.307064               951092890  
131           -0.576034               951097947  
132           -0.462980               951097826  
133           -0.257062               951098075  
134           -0.369061               951097452  

[135 rows x 5 columns]


highly_correlated_neurons = results[results['Highest Correlation'] > 0.7]
highly_correlated_neurons


lowly_correlated_neurons = results[results['Lowest Correlation'] < -0.55]
lowly_correlated_neurons


df = data.copy()
iso_forest = IsolationForest(n_estimators=100, contamination='auto', random_state=42)
iso_forest.fit(df.values)

anomalies = iso_forest.predict(df.values)

df['anomaly_score'] = iso_forest.decision_function(df.values)
df['is_anomaly'] = anomalies

anomalous_neurons = df[df['is_anomaly'] == -1]


neuron_entropies = data.apply(lambda x: entropy(x.value_counts()), axis=1)


# Dataset with anomaly values, entropy, correlation, and neuron with the most correlation
df_entropy_anomaly = pd.DataFrame({
    'Entropy' : neuron_entropies,
    'Anomaly' : df['is_anomaly'],
    'Highest reported correlation with another neuron': highest_correlations.values,
    'Most Correlated Neuron': most_correlated_neurons.values,
    'Variance': df.var(axis=1)
    })
df_entropy_anomaly.sort_values(by='Entropy', ascending=False, inplace=True)
df_entropy_anomaly.reset_index(inplace=True)
type_dict = {'index': int}
df_entropy_anomaly = df_entropy_anomaly.astype(type_dict)


df_entropy_anomaly


df_entropy_anomaly['Anomaly'].value_counts()

 1    116
-1     19
Name: Anomaly, dtype: int64


df_entropy_anomaly.head(20)[df_entropy_anomaly['Anomaly']==-1].shape

/var/folders/_2/p0xm13rx0yq8v664xb_zz8tc0000gn/T/ipykernel_78987/1485505297.py:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  df_entropy_anomaly.head(20)[df_entropy_anomaly['Anomaly']==-1].shape

(14, 6)


df_entropy_anomaly.sort_values(by='Variance', ascending=False).head(20)[df_entropy_anomaly['Anomaly']==-1].shape

/var/folders/_2/p0xm13rx0yq8v664xb_zz8tc0000gn/T/ipykernel_78987/2583719294.py:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  df_entropy_anomaly.sort_values(by='Variance', ascending=False).head(20)[df_entropy_anomaly['Anomaly']==-1].shape

(15, 6)


plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly, x='Entropy', y='Variance', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Entropy vs Variance of Neuron Responses')
plt.xlabel('Entropy')
plt.ylabel('Variance')
plt.legend(title='Anomaly')
plt.show()


df_oracle = pd.read_csv('oracle_df.csv')


df_oracle


def super_oracle(df):
    corr = dict()
    # Total number of images
    total_images = len(df['image_id'].unique())
    # Take a neuron
    for neuron_id in df['neuron_id'].unique():
        responses_list = []
        loom_list = []
        # Image counter
        image_counter = 0
        # Take an image which is shown to the neuron
        for image_id in df[(df['neuron_id'] == neuron_id)]['image_id'].unique():
            # Skip if there is only one trial for the image
            if len(df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id)]['trial_id'].unique()) == 1:
                continue
            # Take a random trial for the image and neuron combination
            random_trial = random.choice(df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id)]['trial_id'].unique())
            # Take the neuron response to the random trial
            response = df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id) & (df['trial_id']==random_trial)]['spike_count'].values[0]
            # Compute the mean over all the other trials for the image
            loom_image = df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id) & (df['trial_id']!=random_trial)]['spike_count'].values.mean()
            # Increment the image counter
            image_counter += 1
            # Append
            responses_list.append(response)
            loom_list.append(loom_image)
        # Compute the correlation between the leave-one-out-mean and the actual response for the neuron on that image
        try:
            correlation = pearsonr(responses_list, loom_list)[0]
        except:
            print(f"Neuron {neuron_id} has a problem")
            return neuron_id
        # Compute the image fraction
        image_fraction = image_counter / total_images
        corr[neuron_id] = correlation*image_fraction
    return pd.DataFrame.from_dict(corr, orient='index', columns=['super_oracle_correlation'])


cor_super = super_oracle(df_oracle)
cor_super


cor_super.reset_index(inplace=True)
cor_super.sort_values(by='super_oracle_correlation', ascending=False, inplace=True)
cor_super


def average_super_oracle(df, iterations=10):
    rank_sums = {}  # Stores sum of ranks
    correlation_sums = {}  # Stores sum of correlations
    total_neuron_ids = df['neuron_id'].unique()

    # Initialize dictionaries
    for neuron_id in total_neuron_ids:
        rank_sums[neuron_id] = 0
        correlation_sums[neuron_id] = 0

    for _ in range(iterations):
        # Run the super_oracle function
        result_df = super_oracle(df)
        if type(result_df) != pd.DataFrame:
            print(f"Error in neuron {result_df}")
            return 0
        # Sort results by correlation in descending order and reset index
        sorted_df = result_df.sort_values(by='super_oracle_correlation', ascending=False).reset_index()

        # Record ranks and sum up correlations
        for index, row in sorted_df.iterrows():
            neuron_id = row['index']
            correlation = row['super_oracle_correlation']
            rank_sums[neuron_id] += index  # index is the rank
            correlation_sums[neuron_id] += correlation

    # Prepare to return the results
    average_ranks = {neuron_id: rank_sums[neuron_id] / iterations for neuron_id in total_neuron_ids}
    average_correlations = {neuron_id: correlation_sums[neuron_id] / iterations for neuron_id in total_neuron_ids}

    # Create DataFrame from the computed averages
    summary_df = pd.DataFrame.from_dict({
        'Average Rank': average_ranks,
        'Average Correlation': average_correlations
    })

    # Sort by average rank
    summary_df = summary_df.sort_values(by='Average Rank')

    return summary_df


av = average_super_oracle(df_oracle, iterations=50)

av


cor_super = av.copy()
cor_super.drop(columns=['Average Rank'], inplace=True)
cor_super.rename(columns={'Average Correlation': 'oracle_correlation'}, inplace=True)
cor_super.reset_index(inplace=True)
cor_super


df_entropy_anomaly_oracle = pd.merge(df_entropy_anomaly, cor_super, left_on='index', right_on='index')
df_entropy_anomaly_oracle.sort_values(by='oracle_correlation', ascending=False, inplace=True)
df_entropy_anomaly_oracle.reset_index(inplace=True)
df_entropy_anomaly_oracle.drop(['level_0'], axis=1, inplace=True)


df_entropy_anomaly_oracle


plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly_oracle[df_entropy_anomaly_oracle['oracle_correlation']>0], x='Entropy', y='oracle_correlation', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Entropy vs Oracle Correlation of Neuron Responses')
plt.xlabel('Entropy')
plt.ylabel('Oracle Correlation')
plt.legend(title='Anomaly')
plt.show()


plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly_oracle, x='Variance', y='oracle_correlation', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Variance vs Oracle Correlation of Neuron Responses')
plt.xlabel('Variance')
plt.ylabel('Oracle Correlation')
plt.legend(title='Anomaly')
plt.show()


df_entropy_anomaly_oracle.sort_values(by='oracle_correlation', ascending= False)


# Take neurons in the 60th percentile of the oracle correlation and 50th percentile of the variance
top_oracle_neurons = df_entropy_anomaly_oracle[(df_entropy_anomaly_oracle['oracle_correlation'] > df_entropy_anomaly_oracle['oracle_correlation'].quantile(0.6)) & (df_entropy_anomaly_oracle['Variance'] > df_entropy_anomaly_oracle['Variance'].quantile(0.5))]
top_oracle_neurons


top_oracle_neurons.to_csv('top_oracle_neurons.csv', index=False)


def laplace():
    return np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]).astype(np.float32)[None, None, ...]

def positive(weight, cache=None):
    weight.data *= weight.data.ge(0).float()
    return cache

class Pyramid(nn.Module):
    _filter_dict = {
        'gauss5x5': np.float32([
            [0.003765, 0.015019, 0.023792, 0.015019, 0.003765],
            [0.015019, 0.059912, 0.094907, 0.059912, 0.015019],
            [0.023792, 0.094907, 0.150342, 0.094907, 0.023792],
            [0.015019, 0.059912, 0.094907, 0.059912, 0.015019],
            [0.003765, 0.015019, 0.023792, 0.015019, 0.003765]]),
        'gauss3x3': np.float32([
            [1 / 16, 1 / 8, 1 / 16],
            [1 / 8, 1 / 4, 1 / 8],
            [1 / 16, 1 / 8, 1 / 16]]
        ),
        'laplace5x5': np.outer(np.float32([1, 4, 6, 4, 1]), np.float32([1, 4, 6, 4, 1])) / 256,

    }

    def __init__(self, scale_n=4, type='gauss5x5', downsample=True, _skip_upsampling=False):
        """
        Setup Laplace image pyramid
        Args:
            scale_n: number of Laplace pyramid layers to construct
            type: type of Gaussian filter used in pyramid construction. Valid options are: 'gauss5x5', 'gauss3x3', and 'laplace5x5'
            downsample: whether to downsample the image in each layer. Defaults to True
        """
        super().__init__()
        self.type = type
        self.downsample = downsample
        self._skip_upsampling = _skip_upsampling
        h = self._filter_dict[type]
        self.register_buffer('filter', torch.from_numpy(h))
        self.scale_n = scale_n
        self._kern = h.shape[0]
        self._pad = self._kern // 2

    def lap_split(self, img):
        N, c, h, w = img.size()
        filter = self.filter.expand(c, 1, self._kern, self._kern).contiguous()

        # the necessary output padding depends on even/odd of the dimension
        output_padding = (h + 1) % 2, (w + 1) % 2

        smooth = F.conv2d(img, filter, padding=self._pad, groups=c)
        if self.downsample:
            lo = smooth[:, :, ::2, ::2]
            if self._skip_upsampling:
                # Technically incorrect implementation of the Laplace pyramid
                lo2 = smooth
            else:
                lo2 = 4 * F.conv_transpose2d(lo, filter, stride=2, padding=self._pad, output_padding=output_padding,
                                         groups=c)
        else:
            lo = lo2 = smooth

        hi = img - lo2

        return lo, hi

    def forward(self, img):
        levels = []
        for i in range(self.scale_n):
            img, hi = self.lap_split(img)
            levels.append(hi)
        levels.append(img)
        return levels

    def __repr__(self):
        return "Pyramid(scale_n={scale_n}, padding={_pad}, downsample={downsample}, type={type})".format(
            **self.__dict__)


class SpatialTransformerPyramid2d(nn.Module):
    def __init__(self, in_shape, outdims, scale_n=4, positive=False, bias=True,
                 init_range=.1, downsample=True, _skip_upsampling=False, type=None):
        super().__init__()
        self.in_shape = in_shape
        c, w, h = in_shape
        #print(f'c: {c}')
        self.outdims = outdims
        self.positive = positive
        self.gauss_pyramid = Pyramid(scale_n=scale_n, downsample=downsample, _skip_upsampling=_skip_upsampling, type=type)
        self.grid = Parameter(torch.Tensor(1, outdims, 1, 2))
        self.features = Parameter(torch.Tensor(1, c * (scale_n + 1), 1, outdims))
        #print(f'self.features first: {self.features.size()}')

        if bias:
            bias = Parameter(torch.Tensor(outdims))
            self.register_parameter('bias', bias)
        else:
            self.register_parameter('bias', None)
        self.init_range = init_range
        self.initialize()

    def initialize(self):
        self.grid.data.uniform_(-self.init_range, self.init_range)
        self.features.data.fill_(1 / self.in_shape[0])

        if self.bias is not None:
            self.bias.data.fill_(0)

    def group_sparsity(self, group_size):
        f = self.features.size(1)
        n = f // group_size
        ret = 0
        for chunk in range(0, f, group_size):
            ret = ret + (self.features[:, chunk:chunk + group_size, ...].pow(2).mean(1) + 1e-12).sqrt().mean() / n
        return ret

    def feature_l1(self, average=True):
        if average:
            return self.features.abs().mean()
        else:
            return self.features.abs().sum()

    def neuron_layer_power(self, x, neuron_id):
        if self.positive:
            positive(self.features)
        self.grid.data = torch.clamp(self.grid.data, -1, 1)
        N, c, w, h = x.size()
        m = self.gauss_pyramid.scale_n + 1
        feat = self.features.view(1, m * c, self.outdims)
        print(f'self.features second: {feat.size()}')
        y = torch.cat(self.gauss_pyramid(x), dim=1)
        y = (y * feat[:, :, neuron_id, None, None]).sum(1)
        return y.pow(2).mean()

    def forward(self, x, shift=None):
        if self.positive:
            positive(self.features)
        self.grid.data = torch.clamp(self.grid.data, -1, 1)
        N, c, w, h = x.size()
        m = self.gauss_pyramid.scale_n + 1
        feat = self.features.view(1, m * c, self.outdims)

        if shift is None:
            grid = self.grid.expand(N, self.outdims, 1, 2)
        else:
            grid = self.grid.expand(N, self.outdims, 1, 2) + shift[:, None, None, :]

        pools = [F.grid_sample(xx, grid, align_corners=True) for xx in self.gauss_pyramid(x)]
        y = torch.cat(pools, dim=1).squeeze(-1)
        y = (y * feat).sum(1).view(N, self.outdims)

        if self.bias is not None:
            y = y + self.bias
        return y

    def __repr__(self):
        c, w, h = self.in_shape
        r = self.__class__.__name__ + \
            ' (' + '{} x {} x {}'.format(c, w, h) + ' -> ' + str(self.outdims) + ')'
        if self.bias is not None:
            r += ' with bias'

        for ch in self.children():
            r += '  -> ' + ch.__repr__() + '\n'
        return r



class Laplace(nn.Module):
    """
    Laplace filter for a stack of data.
    """

    def __init__(self, padding=0):
        super().__init__()
        self._padding = padding
        self.register_buffer('filter', torch.from_numpy(laplace()))

    def forward(self, x):
        return F.conv2d(x, self.filter, padding=self._padding, bias=None)

class LaplaceL2(nn.Module):
    """
    Laplace regularizer for a 2D convolutional layer.
    """

    def __init__(self, padding=0):
        super().__init__()
        self.laplace = Laplace(padding=padding)

    def forward(self, x, weights=None):
        ic, oc, k1, k2 = x.size()
        if weights is None:
            weights = 1.0
        return (self.laplace(x.view(ic * oc, 1, k1, k2)).view(ic, oc, k1, k2).pow(2) * weights).mean() / 2

class Core:
    def initialize(self):
        #log.info('Not initializing anything')
        pass

    def __repr__(self):
        s = super().__repr__()
        s += ' [{} regularizers: '.format(self.__class__.__name__)
        ret = []
        for attr in filter(lambda x: 'gamma' in x or 'skip' in x, dir(self)):
            ret.append('{} = {}'.format(attr, getattr(self, attr)))
        return s + '|'.join(ret) + ']\n'


class Core2d(Core):
    def initialize(self, cuda=False):
        self.apply(self.init_conv)
        if cuda:
            self = self.cuda()

    @staticmethod
    def init_conv(m):
      if type(m) == nn.Conv2d:
          nn.init.xavier_normal_(m.weight.data)
          if m.bias is not None:
             nn.init.constant_(m.bias.data, 0)
    #def init_conv(m):
     #   if isinstance(m, nn.Conv2d):
      #    nn.init.xavier_normal_(m.weight.data)
       #   if m.bias is not None:
        #    m.bias.fill_(0)

# ---------------------- Identity Core ---------------------------
class IdentityCore(nn.Module, Core2d):
    def __init__(self, **kwargs):
        #log.info('Ignoring input {} when creating {}'.format(repr(kwargs), self.__class__.__name__))
        super().__init__()

    def forward(self, x):
        return x

    def regularizer(self):
        return 0.0

# ---------------------- Conv2d Cores with Pooling -----------------------------


class Stacked2dCore_pooling(Core2d, nn.Module):
  def __init__(self, input_channels, hidden_channels, input_kern, hidden_kern, pooling, hidden_pooling, layers=3,
                 gamma_hidden=0, gamma_input=0., skip=0, final_nonlinearity=True, bias=False, skip_nonlin=False,
                 momentum=0.1, pad_input=True, batch_norm=True, laplace_padding=0, laplace_weights_fn=None,**kwargs):
    #log.info('Ignoring input {} when creating {}'.format(repr(kwargs), self._class.name_))
    super().__init__()

    if skip_nonlin:
      print('Skip non-linearity')
    self._input_weights_regularizer = LaplaceL2(padding=laplace_padding)

    self.layers = layers
    self.gamma_input = gamma_input
    self.gamma_hidden = gamma_hidden
    self.input_channels = input_channels
    self.hidden_channels = hidden_channels

    self.skip = skip

    self.features = nn.Sequential()
    self.pooling = pooling
    self.hidden_pooling = hidden_pooling
        # --- first layer
    layer = OrderedDict()
    layer['conv'] = \
            nn.Conv2d(input_channels, hidden_channels, input_kern,
                      padding=input_kern // 2 if pad_input else 0, bias=bias)
    layer['pool'] = nn.AvgPool2d(pooling)
    if batch_norm:
      layer['norm'] = nn.BatchNorm2d(hidden_channels, momentum=momentum)
    if (not skip_nonlin) and (layers > 1 or final_nonlinearity):
      layer['nonlin'] = nn.ELU(inplace=True)
    self.features.add_module('layer0', nn.Sequential(layer))

        # --- other layers
    for l in range(1, self.layers):
       layer = OrderedDict()
       layer['conv'] = \
           nn.Conv2d(hidden_channels if not skip > 1 else min(skip, l) * hidden_channels,
                     hidden_channels, hidden_kern,
                     padding=hidden_kern // 2, bias=bias)
       layer['pool'] = nn.AvgPool2d(hidden_pooling)
       if batch_norm:
           layer['norm'] = nn.BatchNorm2d(hidden_channels, momentum=momentum)
       if (not skip_nonlin) and (final_nonlinearity or l < self.layers - 1):
           layer['nonlin'] = nn.ELU(inplace=True)
       self.features.add_module('layer{}'.format(l), nn.Sequential(layer))
    self.apply(self.init_conv)

    if laplace_weights_fn is not None:
        _, _, h, w = self.features[0].conv.weight.size()
        dist_grid = torch.sqrt(torch.linspace(-1, 1, h)[:, None].pow(2) + torch.linspace(-1, 1, w).pow(2))
        self.register_buffer('laplace_weights', laplace_weights_fn(dist_grid))
    else:
        self.laplace_weights = None


  def forward(self, input_):
    ret = []
    for l, feat in enumerate(self.features):
      do_skip = l >= 1 and self.skip > 1
      input_ = feat(input_ if not do_skip else torch.cat(ret[-min(self.skip, l):], dim=1))
      ret.append(input_)
    return torch.cat(ret, dim=1)
  def laplace(self):
    return self._input_weights_regularizer(self.features[0].conv.weight, weights=self.laplace_weights)
  def group_sparsity(self):
    ret = 0
    for l in range(1, self.layers):
      ret = ret + self.features[l].conv.weight.pow(2).sum(3, keepdim=True).sum(2, keepdim=True).sqrt().mean()
    return ret / ((self.layers - 1) if self.layers > 1 else 1)
  def regularizer(self):
    return self.group_sparsity() * self.gamma_hidden + self.gamma_input * self.laplace()
  @property
  def outchannels(self):
    return len(self.features) * self.hidden_channels




#951092594

class CustomDataset(Dataset):
    def __init__(self, root_dir, dataframe, transform=None, target_size=(60, 96)):
        self.root_dir = root_dir
        self.dataframe = dataframe
        self.transform = transform
        self.target_size = target_size

        # Filter out non-image files and hidden files
        self.image_names = [f for f in os.listdir(self.root_dir) if os.path.isfile(os.path.join(self.root_dir, f)) and not f.startswith('.')]
        self.responses = self.dataframe['951092075'].tolist()

    def __len__(self):
        return len(self.image_names)

    def __getitem__(self, idx):
        image_name = self.image_names[idx]
        image_path = os.path.join(self.root_dir, image_name)

        # Load image
        image = Image.open(image_path).convert('L')  # Convert to grayscale

        # Resize image
        image = Resize(self.target_size)(image)

        # Apply additional transformations
        if self.transform:
            image = self.transform(image)

        # Get corresponding neuronal response
        response = self.responses[idx]

        return image, response






class CustomCNN_2_pooling(nn.Module):
  def __init__(self, input_channels):
    super(CustomCNN_2_pooling, self).__init__()

    # Define the CNN core with three layers
    self.conv_core = Stacked2dCore_pooling(input_channels= input_channels,
                                        hidden_channels=64,
                                        input_kern=4,
                                        hidden_kern=3,
                                        pooling = 2,
                                        hidden_pooling=1,
                                        layers=3,
                                        gamma_hidden=0.1,
                                        gamma_input=0.1,
                                        skip=1,
                                        final_nonlinearity=True,
                                        bias=False,
                                        skip_nonlin=False,
                                        momentum=0.1,
                                        pad_input=True,
                                        batch_norm=True,
                                        laplace_padding=0, stride = 3)

    # Calculate the number of output channels from the CNN core
    self.outchannels = self.conv_core.outchannels
    #self.readout = CustomReadout(self.outchannels)

    self.readout = SpatialTransformerPyramid2d(in_shape=(self.outchannels, 96, 60),
                                                           type = 'gauss5x5',
                                                           scale_n=3,
                                                          # outchannels = 256,
                                                           outdims=1,
                                                           positive=True)
  def forward(self, x):
    # Forward pass through the convolutional core
    features = self.conv_core(x)
    outputs = self.readout(features)
    return outputs


# Here I import a dataset containing all 118 natural scenes of a
mei_dir = '/content/drive/MyDrive/MEI testing/Images_final'
dataframe_mei = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/dataset_for_agne_mean.csv')
dataframe_mei.sort_values(by='frame', inplace=True)
dataframe_mei.set_index('frame', inplace=True, drop=True)
dataframe_mei = dataframe_mei.drop(dataframe_mei.index[0])

# Here I let the images be reduced in 60X96 format
mei_dataset = CustomDataset(mei_dir, dataframe_mei, transform=transforms.ToTensor())


a = np.random.randint(0,117)
plt.imshow(mei_dataset[a][0].squeeze().numpy(), cmap='gray')
plt.axis('off')
plt.show()


neuron_id = '951097947'


# Define the path to the checkpoint file
model_path = '/content/drive/MyDrive/Neuroscience Project/checkpoints_final/checkpoint_epoch_951097947_size_gabor_kernel_59_accuracy.pth'

# Initialize the model
model = CustomCNN_2_pooling(1)  # Replace YourModelClass with the class name of your model

# Load the checkpoint
checkpoint = torch.load(model_path)

# Load the model state dictionary
model.load_state_dict(checkpoint['model_state_dict'])

# Optional: If you have an optimizer saved in the checkpoint and want to load it
# optimizer = YourOptimizerClass(model.parameters(), lr=your_learning_rate)  # Replace YourOptimizerClass with the class name of your optimizer
# optimizer.load_state_dict(checkpoint['optimizer_state_dict'])

# Access other saved information
epoch = checkpoint['epoch']
train_loss = checkpoint['train_loss']
val_loss = checkpoint['val_loss']
train_correlation = checkpoint['train_correlation']
val_correlation = checkpoint['val_correlation']

# Optionally, you can set the model to evaluation mode if needed
model.eval()


# Initialization
c_star_std = 0

for i in range(len(mei_dataset)):
  #std of images
  std = torch.std(mei_dataset[i][0])
  c_star_std += std
# Average over all images
c_star_std /= len(mei_dataset)
# Put in into number format
c_star_std = c_star_std.item()

print(f"Mean standard deviation for our dataset is {c_star_std}")

Mean standard deviation for our dataset is 0.2655758261680603


# Grad computation function
def grad(model, img, i, r, stds):

    # Set the model to evaluation mode (no weights update)
    model.eval()
    # Contrast regularization parameter
    param = 2

    # Ensure image is the correct dimension
    img = img.view(1, 1, 60, 96)
    # Ensure tracking of gradients for the input image
    img.requires_grad_(True)

    # Forward pass through the model
    outputs = model(img)
    output = outputs.sum()

    # Save statistics
    std = torch.std(img)
    r.append(output)
    stds.append(std)

    # Regularize with respect to the computed c_star_std
    output -= param * torch.abs(torch.tensor(c_star_std) - std)

    # Print to see how the optimization is working
    if i % 10 == 0:
      print(f"Loss to minimize is {output} in step {i} with std {std}")

    # Clear previous gradients
    model.zero_grad()

    # Backpropagate to compute gradients w.r.t. input images
    output.backward()

    # Get the gradients of the image and reshape to [H, W]
    image_gradients = img.grad.detach().cpu().squeeze()
    return image_gradients, output #output [60,96] torch tensor



# Define fast fourier transform
def fft_smooth(grad, factor= 1/4):

    if factor == 0:
        return grad.numpy()  # Directly return the numpy array if no smoothing is to be applied

    # Ensure grad is a torch tensor
    if isinstance(grad, np.ndarray):
        grad = torch.from_numpy(grad).float()

    # Obtain dimensions
    h, w = grad.shape

    # Adjust tw for the reduced dimensionality of the rFFT output
    tw = torch.arange(0, w//2 + 1).float()
    th = torch.arange(0, h).float()

    # Compute the frequency weights
    # Distance from origin in frequency space, avoiding division by zero
    radius = torch.sqrt(th[:, None]**2 + tw[None, :]**2) + 1e-8
    t = 1 / radius**factor

    # Normalize the filter to preserve overall energy
    F = t / t.mean()

    # Compute the FFT of the gradient
    grad_fft = torch.fft.rfftn(grad, dim=(0, 1))

    # Apply the filter in the frequency domain
    filtered_fft = grad_fft * F

    # Compute the inverse FFT to transform back to the spatial domain
    filtered_grad = torch.fft.irfftn(filtered_fft, s=(h, w), dim=(0, 1))

    # Convert back to numpy for any subsequent processing
    return filtered_grad.numpy() #[60,96] numpy array



# Scale in a preferred interval [c,d]

def scale(img, c = 0,d = 1):
  # Quantities for scaling
  a = img.min()
  b = img.max()
  # Scaling formula
  sc_matrix =  c + (img - a) * (d - c) / (b - a)
  return np.array(sc_matrix) # [60,96] numpy array


# Generate white Gaussian noise image
# scaled between [0, 1]
height, width = 60,96

wn = scale(np.random.normal(0, 1, size = (height, width)), 0, 1)

#display
plt.imshow(wn, cmap='gray')
plt.axis('off')
plt.show()


# Here we transform the white noise into a tensor [1,1,60,96] readable by our model
scaled_data_tensor = torch.unsqueeze(torch.tensor(normal_data), dim=0).float()
img_tens = torch.unsqueeze(scaled_data_tensor, 0)
print(img_tens.size())

torch.Size([1, 1, 60, 96])


def generate_mei(wn, steps, sigma, dataset= mei_dataset, alpha=0.1, c_star= c_star_std):

  # Here we transform the whote noise into a tensor readable by our model
  scaled_data_tensor = torch.unsqueeze(torch.tensor(wn), dim=0).float()
  img_tens = torch.unsqueeze(scaled_data_tensor, 0)

  # Some statistics:
  # Some Statistics
  r,stds = [],[]

  for i in range(1, steps+1):

    # Compute the gradients for the imput image
    gradients, output = grad(model, img_tens, i, r, stds)
    response = output.item()

    # Transform the input image as a numpy object for smoother computation
    img_tens = img_tens.squeeze().detach().numpy()
    # Fourier Transform to smooth gradients
    fft_grad = fft_smooth(gradients)
    # Do Gradient Ascent step
    img_tens += alpha*fft_grad
    # Perform a Gaussian Blur and reduce sigma at each step
    sigma = sigma + ((0.001 - sigma)*i)/steps
    img_tens = torch.tensor(scipy.ndimage.gaussian_filter(img_tens,sigma,order = 0))
    # Scale to 0-1
    img_tens = scale(img_tens,0,1)
    # Bring it to tensor in order to re-initialize the for loop
    img_tens = torch.unsqueeze(torch.unsqueeze(torch.tensor(img_tens), dim=0).float(), 0)

  final_mei = img_tens.squeeze().numpy()
  #final_response = response
  return final_mei, r, stds


mei_59, r_59, std_59 = generate_mei(wn, steps=1000, sigma=7)

Loss to minimize is 4.375064849853516 in step 10 with std 0.2515137791633606
Loss to minimize is 4.308009624481201 in step 20 with std 0.269641637802124
Loss to minimize is 4.273122310638428 in step 30 with std 0.2788967788219452
Loss to minimize is 4.265032768249512 in step 40 with std 0.2836967408657074
Loss to minimize is 4.283144474029541 in step 50 with std 0.2857617139816284
Loss to minimize is 4.344533920288086 in step 60 with std 0.28649529814720154
Loss to minimize is 4.481871128082275 in step 70 with std 0.2868041396141052
Loss to minimize is 4.77593469619751 in step 80 with std 0.28716763854026794
Loss to minimize is 5.128570079803467 in step 90 with std 0.2877623736858368
Loss to minimize is 5.434421539306641 in step 100 with std 0.2883569896221161
Loss to minimize is 5.6400580406188965 in step 110 with std 0.28283339738845825
Loss to minimize is 5.604672908782959 in step 120 with std 0.26160919666290283
Loss to minimize is 5.5008931159973145 in step 130 with std 0.2428705245256424
Loss to minimize is 5.414088726043701 in step 140 with std 0.2221423089504242
Loss to minimize is 5.378298759460449 in step 150 with std 0.19699597358703613
Loss to minimize is 5.35123348236084 in step 160 with std 0.17513230443000793
Loss to minimize is 5.330989837646484 in step 170 with std 0.15610796213150024
Loss to minimize is 5.315517902374268 in step 180 with std 0.1395726352930069
Loss to minimize is 5.303035736083984 in step 190 with std 0.12524454295635223
Loss to minimize is 5.292255401611328 in step 200 with std 0.11288338899612427
Loss to minimize is 5.282402992248535 in step 210 with std 0.10227282345294952
Loss to minimize is 5.27306604385376 in step 220 with std 0.09321900457143784
Loss to minimize is 5.264244079589844 in step 230 with std 0.08554375916719437
Loss to minimize is 5.256048202514648 in step 240 with std 0.07908265292644501
Loss to minimize is 5.248531341552734 in step 250 with std 0.07368449121713638
Loss to minimize is 5.241790294647217 in step 260 with std 0.06920875608921051
Loss to minimize is 5.235894203186035 in step 270 with std 0.06552410125732422
Loss to minimize is 5.230829238891602 in step 280 with std 0.0625109001994133
Loss to minimize is 5.2265305519104 in step 290 with std 0.060061972588300705
Loss to minimize is 5.222908020019531 in step 300 with std 0.05808281898498535
Loss to minimize is 5.219875812530518 in step 310 with std 0.05649103969335556
Loss to minimize is 5.217343330383301 in step 320 with std 0.05521591380238533
Loss to minimize is 5.215231418609619 in step 330 with std 0.054197508841753006
Loss to minimize is 5.213469505310059 in step 340 with std 0.05338575690984726
Loss to minimize is 5.211986064910889 in step 350 with std 0.05273943394422531
Loss to minimize is 5.21073579788208 in step 360 with std 0.05222510173916817
Loss to minimize is 5.2096781730651855 in step 370 with std 0.0518156997859478
Loss to minimize is 5.20878267288208 in step 380 with std 0.051489558070898056
Loss to minimize is 5.208019256591797 in step 390 with std 0.05122938007116318
Loss to minimize is 5.207365036010742 in step 400 with std 0.051021452993154526
Loss to minimize is 5.206799507141113 in step 410 with std 0.05085493251681328
Loss to minimize is 5.206311225891113 in step 420 with std 0.050721194595098495
Loss to minimize is 5.205882549285889 in step 430 with std 0.05061344429850578
Loss to minimize is 5.205508708953857 in step 440 with std 0.05052636191248894
Loss to minimize is 5.205178260803223 in step 450 with std 0.05045566335320473
Loss to minimize is 5.204885959625244 in step 460 with std 0.05039804428815842
Loss to minimize is 5.204625129699707 in step 470 with std 0.050350870937108994
Loss to minimize is 5.20439338684082 in step 480 with std 0.05031207948923111
Loss to minimize is 5.204188823699951 in step 490 with std 0.05028005689382553
Loss to minimize is 5.204003810882568 in step 500 with std 0.050253450870513916
Loss to minimize is 5.203835964202881 in step 510 with std 0.05023123323917389
Loss to minimize is 5.203686714172363 in step 520 with std 0.05021261051297188
Loss to minimize is 5.203553199768066 in step 530 with std 0.05019690468907356
Loss to minimize is 5.203433513641357 in step 540 with std 0.05018359050154686
Loss to minimize is 5.203327178955078 in step 550 with std 0.05017225444316864
Loss to minimize is 5.203231334686279 in step 560 with std 0.050162557512521744
Loss to minimize is 5.203144550323486 in step 570 with std 0.05015421658754349
Loss to minimize is 5.203067779541016 in step 580 with std 0.05014701932668686
Loss to minimize is 5.202995777130127 in step 590 with std 0.05014077574014664
Loss to minimize is 5.202930927276611 in step 600 with std 0.05013534426689148
Loss to minimize is 5.202873706817627 in step 610 with std 0.05013059824705124
Loss to minimize is 5.202821731567383 in step 620 with std 0.050126466900110245
Loss to minimize is 5.2027740478515625 in step 630 with std 0.05012283846735954
Loss to minimize is 5.202732563018799 in step 640 with std 0.05011964961886406
Loss to minimize is 5.202693939208984 in step 650 with std 0.05011684447526932
Loss to minimize is 5.202657699584961 in step 660 with std 0.05011436715722084
Loss to minimize is 5.202628135681152 in step 670 with std 0.05011218786239624
Loss to minimize is 5.20259952545166 in step 680 with std 0.050110235810279846
Loss to minimize is 5.202574729919434 in step 690 with std 0.05010852962732315
Loss to minimize is 5.202552795410156 in step 700 with std 0.05010700225830078
Loss to minimize is 5.2025322914123535 in step 710 with std 0.05010564625263214
Loss to minimize is 5.202513694763184 in step 720 with std 0.050104446709156036
Loss to minimize is 5.20249605178833 in step 730 with std 0.05010337382555008
Loss to minimize is 5.202479839324951 in step 740 with std 0.050102412700653076
Loss to minimize is 5.2024664878845215 in step 750 with std 0.05010156333446503
Loss to minimize is 5.20245361328125 in step 760 with std 0.050100814551115036
Loss to minimize is 5.2024431228637695 in step 770 with std 0.05010014772415161
Loss to minimize is 5.202430248260498 in step 780 with std 0.05009952932596207
Loss to minimize is 5.202418804168701 in step 790 with std 0.050098977982997894
Loss to minimize is 5.202408313751221 in step 800 with std 0.050098493695259094
Loss to minimize is 5.202400207519531 in step 810 with std 0.050098076462745667
Loss to minimize is 5.202393531799316 in step 820 with std 0.05009772256016731
Loss to minimize is 5.202387809753418 in step 830 with std 0.05009738728404045
Loss to minimize is 5.202383995056152 in step 840 with std 0.05009710043668747
Loss to minimize is 5.202380657196045 in step 850 with std 0.05009683221578598
Loss to minimize is 5.202377796173096 in step 860 with std 0.050096623599529266
Loss to minimize is 5.2023749351501465 in step 870 with std 0.05009642615914345
Loss to minimize is 5.202373027801514 in step 880 with std 0.05009624361991882
Loss to minimize is 5.202371120452881 in step 890 with std 0.050096094608306885
Loss to minimize is 5.202369689941406 in step 900 with std 0.05009595304727554
Loss to minimize is 5.202367782592773 in step 910 with std 0.050095826387405396
Loss to minimize is 5.202367305755615 in step 920 with std 0.050095703452825546
Loss to minimize is 5.202366828918457 in step 930 with std 0.05009560286998749
Loss to minimize is 5.202363967895508 in step 940 with std 0.05009549483656883
Loss to minimize is 5.202361583709717 in step 950 with std 0.05009540915489197
Loss to minimize is 5.2023606300354 in step 960 with std 0.0500953309237957
Loss to minimize is 5.202358722686768 in step 970 with std 0.050095248967409134
Loss to minimize is 5.202357769012451 in step 980 with std 0.050095196813344955
Loss to minimize is 5.202356815338135 in step 990 with std 0.05009514093399048
Loss to minimize is 5.20235538482666 in step 1000 with std 0.0500950887799263


plt.imshow(mei_59, cmap='gray')
plt.axis('off')
plt.show()


# Create a 1x2 grid of subplots
fig, axs = plt.subplots(1, 2, figsize=(10, 5))

# Plot neural response
axs[0].plot(np.array([t.item() for t in r_59]), label='Neural Response')
axs[0].set_title(' Predicted Neural Response')
axs[0].set_xlabel('Steps')
axs[0].set_ylabel('Response Value')

# Plot standard deviation
axs[1].plot(np.array([t.item() for t in std_59]), label='Standard Deviation', color='orange')
axs[1].set_title('Standard Deviation')
axs[1].set_xlabel('Steps')
axs[1].set_ylabel('Standard Deviation')

# Add a legend
#axs[0].legend()
#axs[1].legend()

# Show the plots
plt.tight_layout()
plt.savefig("Statistics.png")
plt.show()


print(f'The maximal neuronal response is {r_59[-1]}')
print(f'The standard deviation of the MEI is {std_59[-1]}')

The maximal neuronal response is 5.203412055969238
The standard deviation of the MEI is 0.05023481696844101


bl = np.zeros((60,96))
plt.imshow(bl, cmap='gray')
plt.axis('off')
plt.show()


mei_check, r_check, std_check = generate_mei(bl, steps=1000, sigma=7)

Predicted regularized response in step 10 is: 3.4271063804626465 with std 0.30462244153022766
Predicted regularized response in step 20 is: 3.278251886367798 with std 0.3353641629219055
Predicted regularized response in step 30 is: 3.264526844024658 with std 0.3388451039791107
Predicted regularized response in step 40 is: 3.2978897094726562 with std 0.3380120098590851
Predicted regularized response in step 50 is: 3.3261239528656006 with std 0.33806490898132324
Predicted regularized response in step 60 is: 3.362445116043091 with std 0.3379369378089905
Predicted regularized response in step 70 is: 3.4305732250213623 with std 0.33763590455055237
Predicted regularized response in step 80 is: 3.5854854583740234 with std 0.3372311294078827
Predicted regularized response in step 90 is: 3.7948713302612305 with std 0.3368215560913086
Predicted regularized response in step 100 is: 4.000244140625 with std 0.3364415764808655
Predicted regularized response in step 110 is: 4.246574401855469 with std 0.3277827501296997
Predicted regularized response in step 120 is: 4.559811115264893 with std 0.31225138902664185
Predicted regularized response in step 130 is: 4.872156620025635 with std 0.2989574074745178
Predicted regularized response in step 140 is: 5.156697750091553 with std 0.28635427355766296
Predicted regularized response in step 150 is: 5.407121181488037 with std 0.27407658100128174
Predicted regularized response in step 160 is: 5.553841590881348 with std 0.2559845447540283
Predicted regularized response in step 170 is: 5.575096607208252 with std 0.23008406162261963
Predicted regularized response in step 180 is: 5.586808204650879 with std 0.20652490854263306
Predicted regularized response in step 190 is: 5.590517997741699 with std 0.1852765530347824
Predicted regularized response in step 200 is: 5.586966037750244 with std 0.16623751819133759
Predicted regularized response in step 210 is: 5.577611923217773 with std 0.14879246056079865
Predicted regularized response in step 220 is: 5.564088344573975 with std 0.1334252953529358
Predicted regularized response in step 230 is: 5.534425735473633 with std 0.11956727504730225
Predicted regularized response in step 240 is: 5.503382205963135 with std 0.10764080286026001
Predicted regularized response in step 250 is: 5.472306728363037 with std 0.0974496528506279
Predicted regularized response in step 260 is: 5.442211151123047 with std 0.0888018012046814
Predicted regularized response in step 270 is: 5.413810729980469 with std 0.08151818066835403
Predicted regularized response in step 280 is: 5.387557029724121 with std 0.07543240487575531
Predicted regularized response in step 290 is: 5.363697528839111 with std 0.07038965076208115
Predicted regularized response in step 300 is: 5.3423261642456055 with std 0.06624595075845718
Predicted regularized response in step 310 is: 5.323410511016846 with std 0.06286844611167908
Predicted regularized response in step 320 is: 5.306825637817383 with std 0.06013665348291397
Predicted regularized response in step 330 is: 5.292392253875732 with std 0.057943157851696014
Predicted regularized response in step 340 is: 5.279891490936279 with std 0.05619363486766815
Predicted regularized response in step 350 is: 5.26910400390625 with std 0.054806943982839584
Predicted regularized response in step 360 is: 5.259816646575928 with std 0.05371422320604324
Predicted regularized response in step 370 is: 5.251826763153076 with std 0.052857816219329834
Predicted regularized response in step 380 is: 5.244958877563477 with std 0.05218995735049248
Predicted regularized response in step 390 is: 5.23906135559082 with std 0.05167163163423538
Predicted regularized response in step 400 is: 5.2339959144592285 with std 0.051271095871925354
Predicted regularized response in step 410 is: 5.229645729064941 with std 0.05096295475959778
Predicted regularized response in step 420 is: 5.225910186767578 with std 0.05072701722383499
Predicted regularized response in step 430 is: 5.222698211669922 with std 0.05054721608757973
Predicted regularized response in step 440 is: 5.219935894012451 with std 0.05041096732020378
Predicted regularized response in step 450 is: 5.217560768127441 with std 0.05030841380357742
Predicted regularized response in step 460 is: 5.215515613555908 with std 0.05023181438446045
Predicted regularized response in step 470 is: 5.213753700256348 with std 0.050175171345472336
Predicted regularized response in step 480 is: 5.212235450744629 with std 0.05013379827141762
Predicted regularized response in step 490 is: 5.21092414855957 with std 0.0501040555536747
Predicted regularized response in step 500 is: 5.209792137145996 with std 0.05008315667510033
Predicted regularized response in step 510 is: 5.208813190460205 with std 0.05006889998912811
Predicted regularized response in step 520 is: 5.2079668045043945 with std 0.05005962401628494
Predicted regularized response in step 530 is: 5.207233905792236 with std 0.050054050981998444
Predicted regularized response in step 540 is: 5.206599235534668 with std 0.05005118250846863
Predicted regularized response in step 550 is: 5.206047058105469 with std 0.05005024001002312
Predicted regularized response in step 560 is: 5.205570697784424 with std 0.05005071312189102
Predicted regularized response in step 570 is: 5.205156326293945 with std 0.05005211755633354
Predicted regularized response in step 580 is: 5.204796314239502 with std 0.0500541515648365
Predicted regularized response in step 590 is: 5.204482078552246 with std 0.050056569278240204
Predicted regularized response in step 600 is: 5.204209804534912 with std 0.050059184432029724
Predicted regularized response in step 610 is: 5.203973293304443 with std 0.0500619076192379
Predicted regularized response in step 620 is: 5.2037672996521 with std 0.05006461217999458
Predicted regularized response in step 630 is: 5.203588485717773 with std 0.050067245960235596
Predicted regularized response in step 640 is: 5.203433036804199 with std 0.050069764256477356
Predicted regularized response in step 650 is: 5.203299045562744 with std 0.05007214471697807
Predicted regularized response in step 660 is: 5.203179359436035 with std 0.05007435008883476
Predicted regularized response in step 670 is: 5.203074932098389 with std 0.050076380372047424
Predicted regularized response in step 680 is: 5.202983856201172 with std 0.05007825419306755
Predicted regularized response in step 690 is: 5.202906131744385 with std 0.05007999762892723
Predicted regularized response in step 700 is: 5.202838897705078 with std 0.05008157715201378
Predicted regularized response in step 710 is: 5.202776908874512 with std 0.050082966685295105
Predicted regularized response in step 720 is: 5.202723026275635 with std 0.05008423700928688
Predicted regularized response in step 730 is: 5.202676773071289 with std 0.05008537694811821
Predicted regularized response in step 740 is: 5.20263671875 with std 0.05008641257882118
Predicted regularized response in step 750 is: 5.202601909637451 with std 0.0500873401761055
Predicted regularized response in step 760 is: 5.20257043838501 with std 0.05008815973997116
Predicted regularized response in step 770 is: 5.202544212341309 with std 0.05008890479803085
Predicted regularized response in step 780 is: 5.202520847320557 with std 0.05008956044912338
Predicted regularized response in step 790 is: 5.202498912811279 with std 0.050090137869119644
Predicted regularized response in step 800 is: 5.202481269836426 with std 0.05009064823389053
Predicted regularized response in step 810 is: 5.202466011047363 with std 0.05009111762046814
Predicted regularized response in step 820 is: 5.202452659606934 with std 0.050091538578271866
Predicted regularized response in step 830 is: 5.202441215515137 with std 0.05009189993143082
Predicted regularized response in step 840 is: 5.2024312019348145 with std 0.05009223520755768
Predicted regularized response in step 850 is: 5.202423572540283 with std 0.05009251832962036
Predicted regularized response in step 860 is: 5.202415943145752 with std 0.05009276419878006
Predicted regularized response in step 870 is: 5.202408790588379 with std 0.050092991441488266
Predicted regularized response in step 880 is: 5.2024030685424805 with std 0.05009319633245468
Predicted regularized response in step 890 is: 5.20239782333374 with std 0.050093356519937515
Predicted regularized response in step 900 is: 5.202394008636475 with std 0.05009353160858154
Predicted regularized response in step 910 is: 5.202390670776367 with std 0.05009365826845169
Predicted regularized response in step 920 is: 5.20238733291626 with std 0.050093770027160645
Predicted regularized response in step 930 is: 5.202384948730469 with std 0.050093889236450195
Predicted regularized response in step 940 is: 5.202383041381836 with std 0.05009398236870766
Predicted regularized response in step 950 is: 5.202381134033203 with std 0.05009406432509422
Predicted regularized response in step 960 is: 5.202378273010254 with std 0.050094123929739
Predicted regularized response in step 970 is: 5.202377796173096 with std 0.05009419843554497
Predicted regularized response in step 980 is: 5.202375888824463 with std 0.05009424313902855
Predicted regularized response in step 990 is: 5.2023749351501465 with std 0.05009428411722183
Predicted regularized response in step 1000 is: 5.20237398147583 with std 0.050094325095415115


plt.imshow(mei_check, cmap='gray')
plt.axis('off')
plt.show()


firing_rate_gabors = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/dataset_for_agne_gabor.csv')

# make it more readable
firing_rate_gabors = firing_rate_gabors.set_index(firing_rate_gabors['stimulus_presentation_id']).drop('stimulus_presentation_id', axis=1)
firing_rate_gabors.head()


firing_rate_gabors[neuron_id].min(), firing_rate_gabors[neuron_id].max()

(0, 8)


gabors_meta = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/gabors_meta.unknown', sep = '\t')
gabors_meta = gabors_meta.set_index('stimulus_presentation_id')   # more readable:
gabors_meta.head()


gabors_meta['duration'].value_counts().index

Index([0.2502017433212131, 0.2502182566787496, 0.2502157566788128,
       0.2502067433212005, 0.2502132566787622, 0.2501967433212257,
       0.2502042433212637, 0.2502207566788001, 0.2502157566786991,
       0.2502182566788065,
       ...
       0.2502017433213837, 0.2502257566789012, 0.2502092433212084,
       0.2502232566786802, 0.2502142433212952, 0.2501992304961362,
       0.2501942433212605, 0.2501917433212668, 0.2502132566787765,
       0.2501917433208973],
      dtype='float64', name='duration', length=230)


print('Sizes:', set(gabors_meta['size']))
print('Phases:', set(gabors_meta['phase']))
print('Spatial Frequencies:', set(gabors_meta['spatial_frequency']))
print('Temporal Frequencies:', set(gabors_meta['temporal_frequency']))
print('Contrast:', set(gabors_meta['contrast']))

Sizes: {'[20.0, 20.0]'}
Phases: {'[3644.93333333, 3644.93333333]'}
Spatial Frequencies: {0.08}
Temporal Frequencies: {4.0}
Contrast: {0.8}


gabors_meta['orientation'].value_counts()

orientation
90.0    1215
0.0     1215
45.0    1215
Name: count, dtype: int64


## CREATE GABORS

def create_gabor(sz, phase, wavelength, orientation, sigma,
                 dy, dx):
    """ Create a gabor patch (sinusoidal + gaussian).

    Arguments:
        height (int): Height of the image in pixels.
        width (int): Width of the image in pixels.
        phase (float): Angle at which to start the sinusoid in degrees.
        wavelength (float): Wavelength of the sinusoid (1 / spatial frequency) in pixels.
        orientation (float): Counterclockwise rotation to apply (0 is horizontal) in
            degrees.
        sigma (float): Sigma of the gaussian mask used in pixels
        dy (float): Amount of translation in y (positive moves down) in pixels/height.
        dx (float): Amount of translation in x (positive moves right) in pixels/height.

    Returns:
        Array of height x width shape with the required gabor.
    """

    # Compute image size to avoid translation or rotation producing black spaces
    padding = max(height, width)
    imheight = height + 2 * padding
    imwidth = width + 2 * padding
    # we could have diff pad sizes per dimension = max(dim_size, sqrt((h/2)^2 + (w/2)^2))
    # but this simplifies the code for just a bit of inefficiency

    # Create sinusoid with right wavelength and phase
    start_sample = phase
    step_size = 360 / wavelength
    samples = start_sample + step_size * np.arange(imheight)
    samples = np.mod(samples, 360)  # in degrees
    rad_samples = samples * (np.pi / 180)  # radians
    sin = np.sin(rad_samples)

    # Create Gabor by stacking the sinusoid along the cols
    gabor = np.tile(sin, (imwidth, 1)).T

    # Rotate around center
    gabor = ndimage.rotate(gabor, orientation, reshape=False)

    # Apply gaussian mask, its is normalized to take values in [0,1]
    gaussy = signal.windows.gaussian(imheight, std=sigma)
    gaussx = signal.windows.gaussian(imwidth, std=sigma)
    mask = np.outer(gaussy, gaussx)
    gabor = gabor * mask

    # Translate (this is only approximate but it should be good enough)
    if abs(dx) > 1 or abs(dy) > 1:
        raise ValueError('Please express translations as factors of the height/width,'
                         'i.e, a number in interval [-1, 1] ')
    dy = int(dy * height)  # int is the approximation
    dx = int(dx * width)
    gabor = gabor[padding - dy: -padding - dy, padding - dx: -padding - dx]

    if gabor.shape != (height, width):
        raise ValueError('Dimensions of gabor do not match desired dimensions.')

    return gabor.astype(np.float32)


def most_exc_gabors(neuron_id, df = firing_rate_gabors):

    neuron_spks = df[neuron_id]

    # Sort the indices based on the firing rates in descending order
    sorted_indices = sorted(neuron_spks.index, key=lambda i: neuron_spks[i], reverse=True)

    # Select the top 3 indices with the highest firing rates
    top_3_indices = sorted_indices[:3]

    return top_3_indices


## Computes the mean gabor as a numpy array

def mean_gabor(neuron_id, gab_meta = gabors_meta):
  gabors_okay = gabors_meta.T.iloc[:,most_exc_gabors(neuron_id)].T
  gabors_okay_ids = list(gabors_okay.index)

  gabors = []
  for i in gabors_okay_ids:
    gabors.append(create_gabor(sz=(60,96), phase=3644.93333333, wavelength=np.pi, orientation=gabors_okay.loc[i,'orientation'], sigma=5,
                   dy=0, dx=0))

  gabors_combined = np.stack(gabors, axis=0)

  # Compute the mean across the arrays along axis 0
  mean_gabors = np.mean(gabors_combined, axis=0)

  # Scale in range [0, 1]
  mean_gabors = scale(mean_gabors, 0, 1)

  return mean_gabors


gabor_59 = mean_gabor(neuron_id)


plt.imshow(gabor_59, cmap='gray')
plt.axis('off')
plt.show()


f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_59, cmap='gray')
axarr[0].axis('off')
axarr[1].imshow(gabor_59, cmap = 'gray');
axarr[1].axis('off')
plt.show()


# Let us check which is the most suitable radius for both gabors and MEIs
x = np.zeros(96)
for i in range(len(x)):
  x[i] = mean_gabors[29][i]

plt.figure(figsize = (20,5))
plt.axvline(x=35, color='red', ls='--')
plt.axvline(x=61, color='red', ls='--')
plt.xticks([i for i in range(1,97)if i%2==0])
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7edd589b3c40>]


x = np.zeros(96)
for i in range(len(x)):
  x[i] = final_mei[29][i]

plt.figure(figsize = (20,5))
plt.axvline(x=40, color='red', ls='--')
plt.axvline(x=66, color='red', ls='--')
plt.xticks([i for i in range(1,97)if i%2==0])
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7edd588d2080>]


def create_circular_mask(h, w, center=None, radius=12):

    if center is None: # use the middle of the image
        center = (int(w/2), int(h/2))
    if radius is None: # use the smallest distance between the center and image walls
        radius = min(center[0], center[1], w-center[0], h-center[1])

    Y, X = np.ogrid[:h, :w]
    dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)

    mask = dist_from_center <= radius
    return mask


masked_59, lum_59 = luminance(mei_59, radius=12, center=(int(mei_59.shape[1]/2+3), int(mei_59.shape[0]/2)))
masked_gab_59, lum_gab_59 = luminance(gabor_59, radius=12, center=None)

f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(masked_59, cmap='gray')
axarr[1].imshow(masked_gab_59, cmap = 'gray');


import math

def luminance(img, radius, center):
  # apply mask:
  h, w = img.shape[:2]
  themask = create_circular_mask(h,w, radius=radius, center=center)
  masked_img = img.copy()
  masked_img[~themask] = 0.0

  # compute luminance on masked img
  lum=0
  for x in masked_img:
    lum += x.sum()
  lum = lum / (math.pi*(radius**2))
  return masked_img, lum


print(f'MEI luminance: {lum_59}')
print(f'Gabor luminance: {lum_gab_59}')

MEI luminance: 0.44963369290239985
Gabor luminance: 0.4097204913589961


print(f'C star: {c_star_std}')

C star: 0.2675706446170807


def masked_contrast(image, mask):
    masked_values = image[mask]       #Apply mask to the image
    std_dev = np.std(masked_values)   #standard deviation
    return std_dev


print("Standard deviation of non-masked MEI:", mei_59.std())
print("Standard deviation of non-masked Gabor:", gabor_59.std())

Standard deviation of non-masked MEI: 0.050091017
Standard deviation of non-masked Gabor: 0.03668658


h, w = 60, 96
radius = 12

mask_mei_c = create_circular_mask(h, w, radius=radius, center= (int(w/2+3), int(h/2)))
masked_contrast_mei = masked_contrast(mei_59, mask_mei_c)
print("Masked Contrast MEI:", masked_contrast_mei)

mask_gab_c = create_circular_mask(h, w, radius=radius)
masked_contrast_gab = masked_contrast(gabor_59, mask_gab_c)
print("Masked Contrast Gabor:", masked_contrast_gab)

Masked Contrast MEI: 0.15114456
Masked Contrast Gabor: 0.13230467


# GABOR
fft_gab = fftshift(fft2(gabor_59))
pow_spec_gab = (np.abs(fft_gab))**2

freq_x = np.fft.fftshift(np.fft.fftfreq(gabor_59.shape[1])) * gabor_59.shape[1]
freq_y = np.fft.fftshift(np.fft.fftfreq(gabor_59.shape[0])) * gabor_59.shape[0]


# MEI
fft_mei = fftshift(fft2(mei_59))
pow_spec_mei = (np.abs(fft_mei))**2

freq_x_mei = np.fft.fftshift(np.fft.fftfreq(mei_59.shape[1])) * mei_59.shape[1]
freq_y_mei = np.fft.fftshift(np.fft.fftfreq(mei_59.shape[0])) * mei_59.shape[0]


## function to compute axial average:
def axial_average(power, bin_size=1/16*np.pi):
    height, width = power.shape
    center = (int(width/2), int(height/2))
    Y, X = np.indices((height, width))

    #polar coordinates
    R = np.sqrt((X - center[0])**2 + (Y - center[1])**2)
    Theta = np.arctan2((Y - center[1]), (X - center[0]))
    theta_bins = np.floor(Theta / bin_size).astype(int)
    theta_bins += -theta_bins.min()  # Shift scale to be positive

    # Average by bins
    radial_sum = np.bincount(theta_bins.ravel(), weights=power.ravel())
    radial_count = np.bincount(theta_bins.ravel())
    radial_mean = radial_sum / radial_count
    return radial_mean, Theta, theta_bins

bin_size = 1/16*np.pi


P_theta, Theta, theta_bins = axial_average(pow_spec_gab)
P_theta_mei, Theta_mei, theta_bins_mei = axial_average(pow_spec_mei)


def sym_index(P_theta, n):
    angles = np.arange(len(P_theta)) * bin_size
    ei2ntheta = np.exp(1j * 2 * n * angles)
    numerator = np.abs(np.sum(P_theta * ei2ntheta))
    denominator = np.sum(P_theta)
    return numerator / denominator


# Calculate for 1-fold and 2-fold symmetry index
## GABOR
SI_1 = sym_index(P_theta, 1)
SI_2 = sym_index(P_theta, 2)
print("Gabor:")
print(f"1-Fold Symmetry Index: {SI_1}")
print(f"2-Fold Symmetry Index: {SI_2}\n")

## MEI
SI_1_mei = sym_index(P_theta_mei, 1)
SI_2_mei = sym_index(P_theta_mei, 2)
print("MEI:")
print(f"1-Fold Symmetry Index: {SI_1_mei}")
print(f"2-Fold Symmetry Index: {SI_2_mei}")

Gabor:
1-Fold Symmetry Index: 0.9992552829452915
2-Fold Symmetry Index: 0.9975716528845747

MEI:
1-Fold Symmetry Index: 0.9696106836897902
2-Fold Symmetry Index: 0.9845564551083268


f, axarr = plt.subplots(1,2, figsize=(30,10))

im0 = axarr[0].imshow(np.log1p(pow_spec_gab), extent=[freq_x.min(), freq_x.max(), freq_y.min(), freq_y.max()], cmap='gnuplot')
f.colorbar(im0, ax=axarr[0])
axarr[0].set_title("Log Power Spectrum of Gabor")
axarr[0].set_xlabel('horizontal freq')
axarr[0].set_ylabel('vertical freq')

im1 = axarr[1].imshow(np.log1p(pow_spec_mei), extent=[freq_x_mei.min(), freq_x_mei.max(), freq_y_mei.min(), freq_y_mei.max()], cmap='gnuplot');
f.colorbar(im1, ax=axarr[1])
axarr[1].set_title("Log Power Spectrum of MEI")
axarr[1].set_xlabel('horizontal freq')
axarr[1].set_ylabel('vertical freq')
plt.show()


neuron_id = '951097919'


f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_46, cmap='gray')
axarr[1].imshow(gabor_46, cmap = 'gray');


neuron_id = '951092949'


f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_34, cmap='gray')
axarr[1].imshow(gabor_34, cmap = 'gray');


neuron_id = '951098928'


f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_32, cmap='gray')
axarr[1].imshow(gabor_32, cmap = 'gray');

	stimulus_presentation_id	unit_id	time_since_stimulus_presentation_onset
spike_time
5909.794764	51355	951097547	0.000318
5909.796721	51355	951093080	0.002274
5909.801697	51355	951098928	0.007251
5909.803021	51355	951092303	0.008574
5909.805297	51355	951098560	0.010851
...	...	...	...
8569.511638	68228	951092303	0.244363
8569.512771	68228	951093688	0.245497
8569.515238	68228	951098021	0.247964
8569.516072	68228	951097682	0.248797
8569.517005	68228	951097869	0.249731

	stimulus_block	start_time	stop_time	frame	stimulus_name	duration	stimulus_condition_id
stimulus_presentation_id
51355	9.0	5909.794447	5910.044666	13.0	natural_scenes	0.250219	4908
51356	9.0	5910.044666	5910.294885	38.0	natural_scenes	0.250219	4909
51357	9.0	5910.294885	5910.545104	30.0	natural_scenes	0.250219	4910
51358	9.0	5910.545104	5910.795324	35.0	natural_scenes	0.250219	4911
51359	9.0	5910.795324	5911.045522	112.0	natural_scenes	0.250198	4912
...	...	...	...	...	...	...	...
68224	13.0	8568.266451	8568.516655	91.0	natural_scenes	0.250204	5005
68225	13.0	8568.516655	8568.766859	83.0	natural_scenes	0.250204	5022
68226	13.0	8568.766859	8569.017064	55.0	natural_scenes	0.250204	5014
68227	13.0	8569.017064	8569.267274	87.0	natural_scenes	0.250211	4955
68228	13.0	8569.267274	8569.517485	105.0	natural_scenes	0.250211	4942

	951092050	951092075	951092303	951092369	951092398	951092410	951092437	951092450	951092475	951092488	...	951098773	951098807	951098850	951098871	951098928	951099272	951099306	951099320	951099475	951099491
frame
-1.0	3.48	2.06	4.54	5.28	0.52	0.68	2.84	1.22	0.16	3.42	...	1.98	2.28	0.08	0.14	4.38	0.10	1.70	0.82	0.18	0.90
0.0	3.14	1.88	2.72	5.46	0.34	2.54	3.54	1.40	0.04	3.38	...	1.62	2.60	0.12	0.04	3.80	0.16	0.80	0.46	0.18	0.74
1.0	3.14	2.18	3.88	5.04	0.72	1.02	4.20	1.04	0.04	2.26	...	1.98	3.84	0.14	0.12	5.08	0.36	0.90	1.30	0.16	2.42
2.0	3.06	1.98	3.38	5.32	3.32	1.68	5.70	1.74	0.04	3.04	...	2.24	3.14	0.10	0.18	6.20	0.14	1.12	0.78	0.08	1.08
3.0	3.16	1.94	3.46	4.34	1.20	2.26	5.28	1.44	0.08	2.48	...	1.50	3.66	0.14	0.08	6.18	0.24	1.00	0.40	0.14	1.28
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
113.0	2.88	1.76	2.40	5.20	2.92	4.12	7.98	1.60	0.04	3.26	...	4.20	5.30	0.30	0.36	6.86	0.10	3.98	1.38	0.02	0.46
114.0	3.14	2.66	3.08	4.64	0.90	0.62	4.62	1.14	0.06	1.98	...	2.18	3.76	0.08	0.18	5.44	0.22	1.32	0.82	0.20	0.80
115.0	3.68	2.02	3.30	5.28	0.68	2.66	6.52	1.36	0.12	2.54	...	2.12	3.70	0.20	0.14	6.00	0.16	2.74	0.64	0.14	1.00
116.0	4.00	2.92	4.00	5.52	0.58	1.30	4.98	1.40	0.12	2.30	...	2.52	2.94	0.10	0.16	4.44	0.14	0.78	0.44	0.06	1.16
117.0	3.54	2.50	3.14	5.16	0.90	1.94	7.12	1.32	0.06	2.04	...	2.68	4.16	0.02	0.10	5.48	0.12	1.28	1.86	0.24	0.98

	frame	stimulus_presentation_id	951092050	951092075	951092303	951092369	951092398	951092410	951092437	951092450	...	951098773	951098807	951098850	951098871	951098928	951099272	951099306	951099320	951099475	951099491
0	-1.0	60469.34	3.48	2.06	4.54	5.28	0.52	0.68	2.84	1.22	...	1.98	2.28	0.08	0.14	4.38	0.10	1.70	0.82	0.18	0.90
1	0.0	58409.40	3.14	1.88	2.72	5.46	0.34	2.54	3.54	1.40	...	1.62	2.60	0.12	0.04	3.80	0.16	0.80	0.46	0.18	0.74
2	1.0	59487.90	3.14	2.18	3.88	5.04	0.72	1.02	4.20	1.04	...	1.98	3.84	0.14	0.12	5.08	0.36	0.90	1.30	0.16	2.42
3	2.0	57120.46	3.06	1.98	3.38	5.32	3.32	1.68	5.70	1.74	...	2.24	3.14	0.10	0.18	6.20	0.14	1.12	0.78	0.08	1.08
4	3.0	59864.54	3.16	1.94	3.46	4.34	1.20	2.26	5.28	1.44	...	1.50	3.66	0.14	0.08	6.18	0.24	1.00	0.40	0.14	1.28
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
114	113.0	59083.06	2.88	1.76	2.40	5.20	2.92	4.12	7.98	1.60	...	4.20	5.30	0.30	0.36	6.86	0.10	3.98	1.38	0.02	0.46
115	114.0	57132.22	3.14	2.66	3.08	4.64	0.90	0.62	4.62	1.14	...	2.18	3.76	0.08	0.18	5.44	0.22	1.32	0.82	0.20	0.80
116	115.0	59239.70	3.68	2.02	3.30	5.28	0.68	2.66	6.52	1.36	...	2.12	3.70	0.20	0.14	6.00	0.16	2.74	0.64	0.14	1.00
117	116.0	58072.78	4.00	2.92	4.00	5.52	0.58	1.30	4.98	1.40	...	2.52	2.94	0.10	0.16	4.44	0.14	0.78	0.44	0.06	1.16
118	117.0	57424.92	3.54	2.50	3.14	5.16	0.90	1.94	7.12	1.32	...	2.68	4.16	0.02	0.10	5.48	0.12	1.28	1.86	0.24	0.98

	stimulus_presentation_id	stimulus_block	start_time	stop_time	stimulus_name	frame	duration	stimulus_condition_id
0	51355	9.0	5909.794447	5910.044666	natural_scenes	13.0	0.250219	4908
1	51356	9.0	5910.044666	5910.294885	natural_scenes	38.0	0.250219	4909
2	51357	9.0	5910.294885	5910.545104	natural_scenes	30.0	0.250219	4910
3	51358	9.0	5910.545104	5910.795324	natural_scenes	35.0	0.250219	4911
4	51359	9.0	5910.795324	5911.045522	natural_scenes	112.0	0.250198	4912
...	...	...	...	...	...	...	...	...
5945	68224	13.0	8568.266451	8568.516655	natural_scenes	91.0	0.250204	5005
5946	68225	13.0	8568.516655	8568.766859	natural_scenes	83.0	0.250204	5022
5947	68226	13.0	8568.766859	8569.017064	natural_scenes	55.0	0.250204	5014
5948	68227	13.0	8569.017064	8569.267274	natural_scenes	87.0	0.250211	4955
5949	68228	13.0	8569.267274	8569.517485	natural_scenes	105.0	0.250211	4942

	index	stimulus_presentation_id	stimulus_block	start_time	stop_time	x_position	phase	stimulus_name	size	y_position	contrast	temporal_frequency	spatial_frequency	orientation	duration	stimulus_condition_id	unit_id	time_since_stimulus_presentation_onset
0	0	1	0.0	84.942787	85.176306	20.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-10.0	0.8	4.0	0.08	90.0	0.233519	1	951093094	0.007928
1	0	1	0.0	84.942787	85.176306	20.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-10.0	0.8	4.0	0.08	90.0	0.233519	1	951098317	0.008801
2	0	1	0.0	84.942787	85.176306	20.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-10.0	0.8	4.0	0.08	90.0	0.233519	1	951092075	0.011728
3	0	1	0.0	84.942787	85.176306	20.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-10.0	0.8	4.0	0.08	90.0	0.233519	1	951092924	0.017895
4	0	1	0.0	84.942787	85.176306	20.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-10.0	0.8	4.0	0.08	90.0	0.233519	1	951098928	0.022368
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
425244	3644	3645	0.0	996.688034	996.938254	-40.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-20.0	0.8	4.0	0.08	0.0	0.250221	123	951092973	0.225382
425245	3644	3645	0.0	996.688034	996.938254	-40.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-20.0	0.8	4.0	0.08	0.0	0.250221	123	951092912	0.225849
425246	3644	3645	0.0	996.688034	996.938254	-40.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-20.0	0.8	4.0	0.08	0.0	0.250221	123	951093795	0.237015
425247	3644	3645	0.0	996.688034	996.938254	-40.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-20.0	0.8	4.0	0.08	0.0	0.250221	123	951093048	0.244382
425248	3644	3645	0.0	996.688034	996.938254	-40.0	[3644.93333333, 3644.93333333]	gabors	[20.0, 20.0]	-20.0	0.8	4.0	0.08	0.0	0.250221	123	951092303	0.249682

unit_id	951092050	951092075	951092303	951092369	951092398	951092410	951092437	951092450	951092475	951092488	...	951098773	951098807	951098850	951098871	951098928	951099272	951099306	951099320	951099475	951099491
stimulus_presentation_id
1	1	4	5	2	1	4	0	1	1	5	...	8	1	1	3	4	0	2	0	0	0
2	1	2	1	2	0	4	0	0	2	1	...	10	2	2	7	5	0	3	0	0	0
3	2	6	2	2	0	8	3	0	1	1	...	5	1	0	6	5	0	2	1	0	0
4	1	3	2	3	1	1	2	1	0	7	...	6	2	0	3	0	0	2	1	0	0
5	3	4	1	3	0	1	15	3	2	0	...	6	4	1	5	3	0	0	0	0	0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3641	3	1	4	3	0	0	0	2	0	0	...	0	0	0	0	0	0	6	0	0	0
3642	3	1	3	2	4	0	0	0	0	0	...	0	0	0	0	2	0	5	0	0	0
3643	3	3	3	2	3	0	1	0	0	0	...	1	0	0	0	2	0	2	0	0	0
3644	2	4	3	8	0	0	0	0	0	0	...	3	1	0	0	1	0	1	0	0	0
3645	2	2	5	2	1	0	0	0	0	0	...	0	0	0	0	0	1	2	2	0	0

	Neuron	Highest Correlation	Most Correlated Neuron	Lowest Correlation	Least Correlated Neuron
32	951093125	0.905601	951098242	-0.553941	951097947
47	951093608	0.736594	951093821	-0.401066	951098279
56	951093821	0.736594	951093608	-0.429395	951092303
60	951093928	0.750732	951098242	-0.526685	951097947
64	951094081	0.848041	951098242	-0.548484	951097947
94	951097826	0.703051	951098560	-0.462980	951099320
105	951098242	0.905601	951093125	-0.654761	951097947
123	951098560	0.703051	951097826	-0.372543	951098303
131	951099306	0.797598	951094081	-0.576034	951097947

	index	Entropy	Anomaly	Highest reported correlation with another neuron	Most Correlated Neuron	Variance
0	951093821	4.465900	-1	0.736594	951093608	2.694286
1	951093608	4.435348	-1	0.736594	951093821	2.418180
2	951092437	4.423117	-1	0.515452	951097547	1.912515
3	951093283	4.400400	-1	0.414326	951092437	1.576702
4	951098928	4.362596	-1	0.503125	951092437	1.434995
...	...	...	...	...	...	...
130	951093810	2.104046	1	0.384974	951094050	0.011959
131	951093510	2.061717	1	0.344498	951094573	0.012383
132	951098317	2.045806	1	0.299732	951093510	0.009211
133	951098372	2.014949	1	0.397437	951093283	0.009025
134	951094036	1.533163	1	0.339248	951099306	0.014631

	super_oracle_correlation
951092050	-0.115015
951092075	-0.016209
951092303	0.134748
951092369	0.089028
951092398	0.517139
...	...
951099272	0.037873
951099475	-0.004608
951098292	0.133259
951094453	0.268118
951098317	-0.029089

	index	super_oracle_correlation
41	951094050	0.651761
108	951093125	0.579137
102	951098242	0.578967
30	951093608	0.569689
4	951092398	0.517139
...	...	...
7	951092488	-0.093438
84	951093416	-0.096325
106	951093020	-0.105021
0	951092050	-0.115015
28	951093565	-0.163587

	Average Rank	Average Correlation
951098242	1.44	0.602705
951093125	2.94	0.571175
951092798	5.70	0.489572
951093846	6.88	0.499745
951093821	7.84	0.452857
...	...	...
951098317	116.82	-0.043106
951093858	116.84	-0.044025
951093565	117.58	-0.081306
951093416	122.58	-0.082603
951092673	126.20	-0.085904

	index	Entropy	Anomaly	Highest reported correlation with another neuron	Most Correlated Neuron	Variance	oracle_correlation
0	951098242	3.992573	1	0.905601	951093125	0.968955	0.602705
1	951093125	3.716815	1	0.905601	951098242	0.497883	0.571175
2	951093846	3.899434	1	0.677264	951092450	0.769341	0.499745
3	951092798	4.129662	1	0.482546	951092949	0.725587	0.489572
4	951093821	4.465900	-1	0.736594	951093608	2.694286	0.452857
...	...	...	...	...	...	...	...
130	951093858	2.302911	1	0.368254	951094654	0.011904	-0.044025
131	951093565	3.004287	1	0.279689	951094050	0.016774	-0.081306
132	951093416	3.134690	1	0.436531	951099306	0.027655	-0.082603
133	951092673	2.696601	1	0.410253	951097530	0.019744	-0.085904
134	951098372	2.014949	1	0.397437	951093283	0.009025	NaN

	index	Entropy	Anomaly	Highest reported correlation with another neuron	Most Correlated Neuron	Variance	oracle_correlation	Color
0	951098242	3.992573	1	0.905601	951093125	0.968955	0.602705	Other Neurons
1	951093125	3.716815	1	0.905601	951098242	0.497883	0.571175	Other Neurons
2	951093846	3.899434	1	0.677264	951092450	0.769341	0.499745	Other Neurons
3	951092798	4.129662	1	0.482546	951092949	0.725587	0.489572	Other Neurons
4	951093821	4.465900	-1	0.736594	951093608	2.694286	0.452857	Other Neurons
5	951097947	4.232909	-1	0.368022	951093608	0.759370	0.446460	Marti Neuron
6	951094050	4.207487	-1	0.529467	951092695	0.848769	0.435930	Other Neurons
7	951092398	3.808470	1	0.655551	951092410	0.378662	0.428909	Other Neurons
8	951093928	4.130337	1	0.750732	951098242	0.877427	0.418735	Other Neurons
9	951094081	3.507181	1	0.848041	951098242	0.603391	0.403576	Other Neurons
10	951093608	4.435348	-1	0.736594	951093821	2.418180	0.400721	Other Neurons
11	951093725	2.987012	1	0.675827	951092940	0.228417	0.364932	Other Neurons
12	951093662	3.655037	1	0.602763	951094559	0.284722	0.353167	Other Neurons
13	951099306	3.979247	1	0.797598	951094081	0.391902	0.348804	Other Neurons
14	951092949	4.111789	-1	0.611074	951093608	1.177489	0.323257	Other Neurons
15	951092940	4.349405	-1	0.675827	951093725	2.236061	0.320312	Other Neurons
16	951093283	4.400400	-1	0.414326	951092437	1.576702	0.313524	Other Neurons
17	951094516	4.148279	1	0.474654	951092949	0.410199	0.295350	Other Neurons
18	951092973	4.301265	-1	0.559254	951097547	1.121742	0.295022	Other Neurons
19	951092410	4.041892	1	0.655551	951092398	0.579999	0.293418	Other Neurons
20	951093349	3.616517	1	0.482599	951092410	0.193406	0.291107	Other Neurons
21	951098928	4.362596	-1	0.503125	951092437	1.434995	0.286298	Other Neurons
22	951092924	3.903790	1	0.252571	951098384	0.158347	0.279925	Other Neurons
23	951093472	4.131069	1	0.497556	951093662	0.828043	0.279441	Other Neurons
24	951092961	4.271673	-1	0.478252	951092983	0.805439	0.276420	Other Neurons
25	951093688	4.336214	1	0.672685	951093608	0.824899	0.271322	Other Neurons
26	951097826	4.057939	1	0.703051	951098560	0.241768	0.260815	Other Neurons
27	951098332	3.943135	1	0.363750	951097653	0.279123	0.259099	Other Neurons
29	951097708	3.697292	1	0.559016	951093846	0.174054	0.239430	Other Neurons
30	951092450	3.951929	1	0.677264	951093846	0.244408	0.235229	Other Neurons
31	951092983	4.170997	1	0.539202	951097919	0.472773	0.228897	Other Neurons
32	951093484	3.909350	1	0.397126	951094529	0.170366	0.224352	Other Neurons
33	951097504	4.165286	1	0.565490	951092450	0.556476	0.223981	Other Neurons
34	951097919	4.069588	1	0.539202	951092983	0.254956	0.217594	Other Neurons
35	951097653	3.692260	1	0.617355	951098360	0.231239	0.213991	Other Neurons
36	951097934	3.846258	1	0.691230	951093125	0.196590	0.204507	Other Neurons
37	951099320	3.872860	1	0.462756	951098242	0.191630	0.204198	Other Neurons
38	951098105	4.005780	1	0.485009	951094081	0.360385	0.202622	Other Neurons
39	951098075	3.906645	1	0.380846	951093928	0.310396	0.201916	Other Neurons
40	951097857	3.859465	1	0.582236	951093662	0.156676	0.193375	Other Neurons

Investigating Optimal Visual Input for Cortical Neurons:¶

A CNN-based MEI approach compared to Gabor patches¶

Authors:¶

INTRODUCTION¶

SCOPE OF THE STUDY:¶

METHODOLOGY:¶

Imports:¶

Creation of the dataset¶

Dataset for CNN and MEI:¶

Dataset for Oracle Correlation Analysis:¶

Most Important Neurons (MINs)¶

Rasterplot analysis¶

Covariance analysis:¶

Isolation Forest for outlier detection¶

Entropy and variance plots:¶

Oracle Correlation¶

CNN Model¶

CNN Structure¶

Core block¶

Readout Layer¶

Most Excitatory Inputs (MEIs):¶

Setup: Dataset and Model¶

Target contrast:¶

MEI computation: The Algorithm¶

Utils:¶

1. Gaussian white noise¶

2.- 5. Iterative process:¶

Statistics Analysis:¶

Increasing neuronal response:¶

Decreasing Standard Deviation:¶

Further Considerations on Standard Deviation:¶

Invariance with respect to initial input:¶

Gabor Filters¶

Preliminary discussion:¶

Generate most responsive gabors:¶

Preliminary analysis:¶

Gabor extraction:¶

Comparative Analysis:¶

Validating MEIs¶

Initial Observations:¶

Orientation:¶

Spatial Frequency:¶

Metrics:¶

Mask:¶

Luminance¶

Contrast:¶

Symmetry:¶

Log Power Spectrum and Spatial Frequency Content:¶

Conclusions:¶

Limitations and Further Analysis:¶

APPENDIX: Other Sampled Neurons¶

NEURON 951097919, accuracy value: $46$%¶

NEURON 951092949, accuracy value: $34$%¶

NEURON 951098928, accuracy value: $32$%¶